clear

theta = [0.5 : 1. : 180.]';
phi   = [0.5 : 1. : 360.]';

File_1 = 'C:\Manoj\ocean\manoj_2006\1x1_summer_2002_new.res';
%File_2 = 'apr00_run_2_new_PC.res';

%File_1 = '200004_run_1.res';
%File_2 = '200004_run_2.res';

fid = fopen(File_1, 'r', 'n');

LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double'); 
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock, 'bof'); % LBlock if file is written on sun

LBlock
np
nt

Hpr = fread(fid, [nt np], 'float');
Hpi = fread(fid, [nt np], 'float');
Htr = fread(fid, [nt np], 'float');
Hti = fread(fid, [nt np], 'float');
Hrr = fread(fid, [nt np], 'float');
Hri = fread(fid, [nt np], 'float');

Epr = fread(fid, [nt np], 'float');
Epi = fread(fid, [nt np], 'float');
Etr = fread(fid, [nt np], 'float');
Eti = fread(fid, [nt np], 'float');

Hx1  = -(Htr+j*Hti)*100*4.*pi;
Hy1  =  (Hpr+j*Hpi)*100*4.*pi;
Hr1  =  (Hrr+j*Hri)*100*4.*pi;

Ex1  = -(Etr+j*Eti)*1.E+06;
Ey1  =  (Epr+j*Epi)*1.E+06;

fclose(fid);

fid = fopen(File_2, 'r', 'n');

LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double'); 
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock, 'bof'); % LBlock if file is written on sun

Hpr = fread(fid, [nt np], 'float');
Hpi = fread(fid, [nt np], 'float');
Htr = fread(fid, [nt np], 'float');
Hti = fread(fid, [nt np], 'float');
Hrr = fread(fid, [nt np], 'float');
Hri = fread(fid, [nt np], 'float');

Epr = fread(fid, [nt np], 'float');
Epi = fread(fid, [nt np], 'float');
Etr = fread(fid, [nt np], 'float');
Eti = fread(fid, [nt np], 'float');

Hx2  = -(Htr+j*Hti)*100*4.*pi;
Hy2  =  (Hpr+j*Hpi)*100*4.*pi;
Hr2  =  (Hrr+j*Hri)*100*4.*pi;

Ex2  = -(Etr+j*Eti)*1.E+06;
Ey2  =  (Epr+j*Epi)*1.E+06;

fclose(fid);

map = [    0.2510         0    0.2510
    0.2331         0    0.3045
    0.2151         0    0.3580
    0.1972         0    0.4115
    0.1793         0    0.4650
    0.1613         0    0.5185
    0.1434         0    0.5720
    0.1255         0    0.6255
    0.1076         0    0.6790
    0.0896         0    0.7325
    0.0717         0    0.7860
    0.0538         0    0.8395
    0.0359         0    0.8930
    0.0179         0    0.9465
         0         0    1.0000
    0.0091    0.0909    1.0000
    0.0182    0.1818    1.0000
    0.0273    0.2727    1.0000
    0.0364    0.3636    1.0000
    0.0455    0.4545    1.0000
    0.0545    0.5455    1.0000
    0.0636    0.6364    1.0000
    0.0727    0.7273    1.0000
    0.0818    0.8182    1.0000
    0.0909    0.9091    1.0000
    0.1000    1.0000    1.0000
    0.2500    1.0000    1.0000
    0.4000    1.0000    1.0000
    0.5500    1.0000    1.0000
    0.7000    1.0000    1.0000
    0.8500    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    0.8333
    1.0000    1.0000    0.6667
    1.0000    1.0000    0.5000
    1.0000    1.0000    0.3333
    1.0000    1.0000    0.1667
    1.0000    1.0000         0
    1.0000    0.9333         0
    1.0000    0.8667         0
    1.0000    0.8000         0
    1.0000    0.7333         0
    1.0000    0.6667         0
    1.0000    0.6000         0
    1.0000    0.5333         0
    1.0000    0.4667         0
    1.0000    0.4000         0
    1.0000    0.3333         0
    1.0000    0.2667         0
    1.0000    0.2000         0
    1.0000    0.1333         0
    1.0000    0.0667         0
    1.0000         0         0
    1.0000         0    0.0909
    1.0000         0    0.1818
    1.0000         0    0.2727
    1.0000         0    0.3636
    1.0000         0    0.4545
    1.0000         0    0.5455
    1.0000         0    0.6364
    1.0000         0    0.7273
    1.0000         0    0.8182
    1.0000         0    0.9091
    1.0000         0    1.0000];
colormap(map);

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 20.5 29]); 

load coast;

subplot(2,1,1);
%imagesc(phi, theta, real(Hr1), [-4 4] );
imagesc(phi, theta, real(Hy1), [-4000 4000] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k'); hold off
colorbar('h');
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
%title 'Ex; mV/km; run = 1; apr00';
title 'Hr; nT; run = 1; apr00';

subplot(2,1,2);
%imagesc(phi, theta, real(Hr2), [-4 4] );
imagesc(phi, theta, real(Hy2), [-4000 4000] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k'); hold off
colorbar('h');
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
%title 'Ex; mV/km; run = 2; apr00';
title 'Hr; nT; run = 2; apr00';

%print -depsc -tiff 'Re_Hr_r1_r2_apr00.eps';

return

fid = fopen('Ex_Ey_3000.txt', 'w');
fprintf(fid, '  M2 tide; rho_lithosphere = 3000; x - north, y - east  \n');
fprintf(fid, '\n');
fprintf(fid, ' theta, phi, ReEx, ImEx, ReEy, ImEy \n');
for iy = 1 : nt
    for ix = 1 : np
        fprintf(fid, '%14.5e    ',  theta(iy));
        fprintf(fid, '%14.5e    ',    phi(ix));
        fprintf(fid, '%14.5e    ',  real(Ex1(iy,ix)));
        fprintf(fid, '%14.5e    ',  imag(Ex1(iy,ix)));
        fprintf(fid, '%14.5e    ',  real(Ey1(iy,ix)));
        fprintf(fid, '%14.5e    ',  imag(Ey1(iy,ix)));
        fprintf(fid, '\n');
    end
end
fclose(fid);

fid = fopen('Ex_Ey_300.txt', 'w');
fprintf(fid, '  M2 tide; rho_lithosphere = 300; x - north, y - east  \n');
fprintf(fid, '\n');
fprintf(fid, ' theta, phi, ReEx, ImEx, ReEy, ImEy \n');
for iy = 1 : nt
    for ix = 1 : np
        fprintf(fid, '%14.5e    ',  theta(iy));
        fprintf(fid, '%14.5e    ',    phi(ix));
        fprintf(fid, '%14.5e    ',  real(Ex2(iy,ix)));
        fprintf(fid, '%14.5e    ',  imag(Ex2(iy,ix)));
        fprintf(fid, '%14.5e    ',  real(Ey2(iy,ix)));
        fprintf(fid, '%14.5e    ',  imag(Ey2(iy,ix)));
        fprintf(fid, '\n');
    end
end
fclose(fid);
